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Abstract. We develop an effective numerical method of studying large-time properties of reversible 
reaction-diffusion systems of type A -f B ^ C with initially separated reactants. Using it we find that there 
are three types of asymptotic reaction zones. In particular we show that the reaction rate can be locally 
negative and concentrations of species A and B can be nonmonotonic functions of the space coordinate x, 
locally significantly exceeding their initial values. 
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1 Introduction 

■ Dynamic reaction fronts formed between initially sepa- 
, rated reactants A and B that perform Brownian motion 
' and react upon contact are an important component of 
\ many physical, chemical and biological systems QIH]. Most 

■ and experimental I16''17''18",'I9'"2U' research has been fo- 
, cused on irreversible reactions of type A -I- B — > C, which 
' exhibit many unexpected phenomena. For example, the 
' width of the reaction zone grows with time t as with 
> surprisingly small value of a = 1/6 [HE], the center of the 

■ reaction front can spontaneously change the direction of 
, its motion [^17| , and the mean-field approximation of the 

■ local reaction rate breaks down at and below the critical 
I dimension dc = 2 [51l7ltT3] . 

In reality, however, most chemical reactions are rever- 
' sible. The simplest model of such a system |23 is based on 
. an assumption that concentrations a, h, and c of species 
' A, B, and C, respectively, effectively depend on time t 
' and only one space coordinate x (even though the system 
. is three-dimensional), and their evolution is governed by 
' three reaction-diffusion equations 

I 

da{x,t) d^a{x,t) 

= (1) 

db{x,t) d%(x,t) , , 

dc(x,t) „ d^c(x,t) , , 
where the effective local reaction rate R{x, t) equals to the 

k 

difference between the production (A + B ^ C) and decay 
^ e-mail: zkozaOift . uni . wroc. pi 



(C ^ A + B ) rates of species C, 

R{x, t) = ka{x, t)b{x, t) — gc{x, t). (4) 

Here Da, Db, and Dc are diffusion coefficients of species 

A, B, and C, respectively, and k,g > are reaction rate 
constants. It is also assumed that initially species A and B 
are uniformly distributed on opposite sides of x = with 
concentrations oq and bo, respectively, 

a{x, 0) = aoH{x), b{x, 0) = boH{-x), c{x, 0) = 0,(5) 

where H{x) is the Heaviside step function (which is for 
X < and 1 for x > 0). Such an initial condition is of- 
ten adopted in experiments jlfilll7lll8lll9ll2()| and simpli- 
fies theoretical analysis, as it enables reduction of a three- 
dimensional problem to a one-dimensional one. 

This model was first studied by Chopard et al. pT] . 
They found that (a) the front width of a reversible reaction 
asymptotically scales with time as if the process was gov- 
erned solely by diffusion {w{t) cx t^/^) and (b) the mean- 
field approximation can be safely applied for systems 
of spatial dimension d — 1,2, 3. However, the fundamental 
problem of giving a detailed description of spatiotemporal 
evolution of reversible reaction-diffusion systems remained 
open until quite recently. 

This problem was recently considered by Sinder and 
Pelleg |11II12II13| . They focused their attention mainly on 
the limit of a vanishingly small backward reaction rate g 
and found that in this limit concentrations of species A, 

B, and C assume the forms typical of irreversible reac- 
tions {g — Q) everywhere except in a very narrow reaction 
zone. They confirmed the result of Ref. that there is 
a crossover between intermediate-time "irreversible" and 
large-time "reversible" regimes. They showed that the asymp- 
totic reaction rate R can have one or two maxima and 
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can even be locally negative (for irreversible reactions R 
always has a single maximum an can never be negative). 
Moreover, they presented strong arguments supporting a 
conjecture that reversible reaction-diffusion processes can 
be divided into two distinct universality classes. One of 
them contains systems with immobile reaction product C 
and asymptotically immobile reaction front, while systems 
with all other combination of the control parameters form 
the other universality class. 

In our recent paper [221 developed a new approach, 
enabling one to analyze the large-time limit of reversible 
reaction-diffusion systems directly, without having to solve 
the original partial differential equations Q - © and 
then taking the limit t — > oo. We proved that in the large- 
time limit functions a{x,t), b{x,t), c{x,t), and R{x,t) ef- 
fectively depend on x only through ^ = xj \/t and take on 
a form 

a[x,t) = A{5,), h{x,i) = B{i), c{x,t)=CiO, 

R{x,t)^t-^ni{0 (6) 

where the scaling functions A, B, C, and TZi are completely 
determined by four equations 



Da 



de 



kAB - gC = 0, 
Til 



1 dA 
2^~d^ 



1 dB 



^ d^C IdC 
with the boundary conditions 



7^l 



lim AiO = flo, lini A{0 = 0, 
lim 6(0 = 0, lim 6(0 = bo, 



(7) 

(8) 

(9) 
(10) 

(11) 
(12) 



Compared with the original problem of solving Eqs. ^ 
- ©, this new approach has two advantages. First, it in- 
volves only ordinary differential equations. Second, it per- 
tains directly to the large-time limit. 

In principle equations - (|10|l completely determine 
the asymptotic, large-time spatiotemporal evolution of an 
arbitrary reversible reaction-diffusion system. Unfortunately, 
they are quite complex and a complete analytical solution 
is known only for the case Da = Db = Dc |221- The aim 
of our paper is to examine these equations numerically for 
other values of the control parameters. 



2 Numerical results 

By measuring length, time, and concentration in units of 
y/DA/kao, 1/fcao, and uq, respectively, the general prob- 
lem of solving - l(Tn|l for arbitrary values of ag, bo. Da, 
Db, Dc, k, and g can be reduced to the one with fIT\ 



We shall adopt these particular values in our further anal- 
ysis. This will leave us with four independent control pa- 
rameters: g, bo, Db, and Dc- 

Our basic equations Q - l(Tn|l can be reduced to two 
ordinary differential equations with two unknown func- 
tions AiO and 6(0, 

d^{A + Dcg-'AB) ^ l^ d{A + g-'AB) 

de 2^ d^ ^ ' 

d\DBB + Dcg-^AB) _ I d{B + g-^AB) 
le ~"2^ dl 

To solve them we employed an iterative method. We first 
assumed that Bo{C) — and, using standard techniques 
'23J , solved H14(l as a linear ordinary differential equation 
for Ao (0 with boundary condition Hll|) . We inserted this 
solution into H15() , which was then solved as a linear differ- 
ential equation for Bo{C)- This solution was again inserted 
into 114|) and used to determine the next approximation 
of Ao{£,)- This procedure was repeated until a required 
accuracy was achieved. 

Taking 6o(0 = as the first approximation leads to 
exact solution for =0 (or, equivalently. A; = 0) after 
the first iteration cycle, and ensures quick convergence for 
most choices of system parameters except when g 1. 
In this case the reaction zone is very narrow and inside it 
functions Ao (0 and Bq (0 vary very rapidly, which makes 
the direct iterative method unstable. This problem may 
be circumvented by first solving (|14|1 and (|15|l for g ^ 1 
and then decreasing g gradually until it reaches the re- 
quired value, each time employing a solution obtained for 
larger g as an initial guess for a smaller value of g. Once 
.4(0 and 6(0 have been determined, the two remaining 
functions of primary interest, C(0 and 7?,i(0 can be cal- 
culated directly from (TJl and ©, respectively. 

The iterative method cannot be applied directly for 
Db — Dc ~ 0, as in this case the left-hand side of p5|l 
vanishes and the order of this differential equation equals 
1 rather than 2. Nevertheless, since in this very particular 
case b{x,t) + c{x,t) = boH{x), after some simple algebra 
we can reduce (|14|) and (|15|l to a single equation 



<PA{i) 

dO 



I dM^ 

2^ 



feog(0 



(16) 



Da - 1, 



(13) 



Although this equation looks very complicated, it can be 
solved quite easily through standard numerical methods. 

To estimate accuracy of the iterative method, we used 
it to solve equations (|14|) and H15|) for the case of equal 
diffusion constants. Da = Db = Dc — 1, and compared 
the results with the exact solutions obtained in Ref. 221 ■ 
For bo = 1, 0.1, 0.01, g = 100, 1, 0.01, and -5 < C < 5 we 
found the relative error to be less than 10^^. 

Next we used (|14|l - (|16f) to investigate thoroughly var- 
ious combinations of system parameters. To ensure that 
the boundary conditions (|ll|l - H12(l are actually satisfied, 
we used a rather wide range —10 < ^ < 10 and com- 
pared the results thus obtained with those calculated for 
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Fig. 1. A, B, C, and TZ\ (asymptotic concentrations of species 
A, B, C, and the scaling function of the reaction rate, re- 
spectively), as functions of ^ = x/\/t for ao = 1, 60 = 0.5, 
Da = Db = Dc = 1, fc = 1 and g = 0.01. Arbitrary units. 

— 15 < C ^ 15, finding no significant differences. More- 
over, upon a thorough numerical scanning of the four- 
dimensional parameter space we found that the solutions 
of equations Q - H1U|I are continuous functions of Db, 
Dc, g, and bo (even when going from one of Sinder and 
Pelleg's universality classes to the other), and can be di- 
vided into three major categories distinguished by specific 
forms of the local reaction rate TZi . 



2.1 Reaction fronts of type I 

A characteristic feature of reaction fronts of type I is that 
the asymptotic reaction rate Tli{£,) is positive for all ^ 
and has a single maximum, which may be identified with 
the reaction front center C/- A typical example of such a 
reaction front is illustrated in Fig. ^ which was obtained 
for bo = 0.5, Da = Db = Dc ^ 1, and g = 0.01. As 
in this case ^/ > 0, we may say that the reaction front 
moves towards the right-hand side of the system. This 
type of solution always appears for Da ^ Db — Dc [221 
and for 5 = j3j, and so we expect it also to appear for 
Da^ Db^^ Dc or for g < 1. 

Interestingly, it turns out that the reaction front formed 
in the case Db = Dc = also belongs to this category. 
This is clearly seen in Fig. [21 obtained for = Dc — 0, 
bo = 0.5, and g = 0.02. A characteristic feature of this 
case is discontinuity of B{£,) and C(^) at f = 0. This re- 
flects the presence of the Heaviside function in Eq. 



2.2 Reaction fronts of type II 

An example of the second type of the asymptotic solution 
is depicted in Fig.[31 which was obtained for much smaller 
value of Dc = 0.01 (the values of other control parameters 
were Db = 0.5, bo ~ 0.25, and g = 0.01). In this case 
7^l has two maxima £_'^''''' w -0.26 and Ci"*"" ~ 1-38. As 
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Fig. 2. A, B, C, and TZi as functions of ^ = xjyft for ao = 1, 
60 = 0.5, Da = 1, -Di3 = I>c = 0, fc = 1 and g = 0.02. Note 
that Z?(^) and C(^) vanish for ^ < 0. Arbitrary units. 
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Fig. 3. A; B, C, and TZi as functions of ^ = xj^ft for ao = 1, 
60 = 0.25, Da = Db = 1, Dc = 0.01, fc = 1 and g = 0.01. 
Note that 7li{^) can assume negative values. Arbitrary units. 

they are of opposite signs, the system apparently has two 
reaction fronts moving at opposite directions. Moreover, 
the reaction rate has one minimum, « 0.19, at which 
it attains a negative value. In the region where TZi < the 
backward reaction (C A -I- B) is thus locally faster than 

the forward reaction {A + B C), although of course the 
global reaction rate T^iiO > 0. 

Formation of a region with negative value of TZi can be 
understood as follows. Consider an asymmetric reaction- 
diffusion system with very small diffusion coefficient of 
species C {Dc ^ Da,Db) and a small backward reac- 
tion rate g. As the reaction proceeds, the reaction front 
moves through the system, leaving behind a region filled 
with practically immobile and very slowly decaying re- 
action product C. At some moment the mobile reaction 
front will leave this region, and so the backward reaction, 
however small, may start to dominate. This may lead to 
formation of a region where TZii^) attains a locally mini- 
mal, perhaps even negative value. 

The dominant backward reaction should lead to pro- 
duction of additional molecules of type A and B. Because 
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Fig. 4. A, B, C, and TZi as functions of ^ = xj^ft for ao — 1, 
feo = 0.25, Da = 1, Db ^ 32, i3c = 0.01, A; = 1 and g = 0.01. 
Arbitrary units. 
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Fig. 5. A, B, C, and TZi as functions of ^ = xjyft for ao = 1, 
feo = 0.5, Da = 1, Db = 0.1, Dc = 10, fc = 1 and g = 1. 
Note that .4o(C) and So (5) are nonmonotonic and 72.i can be 
negative. Arbitrary units. 



Da ,Db ^ Dc , these molecules can easily diffuse away 
from the region filled with molecules C. Then, at the other 
edge of the region rich in species C, they should give a sig- 
nificant contribution to the forward reaction rate, forcing 
T^iiO to change its sign back to positive and forming the 
other reaction front. Such scenario is confirmed by Fig. 13 
which shows that molecules B are present in the whole re- 
gion densely occupied by molecules C, including a region 
between ^^'^^ and 0. Molecules of type B are present in 
this region even though the main reaction front, located 
near is constantly moving away. 

Reaction fronts of type II were first observed by Sinder 
and Pelleg "TT in systems with Dc = 0. They came to the 
conclusion that for mobile reaction fronts {xf{t) ^ 0) the 
larger maximum is located near the point where a(x, t) k, 
b{x,t). However, we found that the opposite situation is 
also possible. This is illustrated in Fig. 01 obtained for 
Dc = 0.01, Db = 32, bo = 0.25, and g = 0.01. As we 
can see, in this case a{x) « b{x) near the second, much 
smaller maximum of Ri (x) . 

We expect reaction fronts of type II to be typical of 
systems where either Db or Dc are much smaller than 
Da- We base this conjecture on Eqs ((HJ and pOfl which 
ensure that if Dc = or Db = then i?i(0) = 0. Since 
i?i(^) is continuous we may thus expect that at least for 
highly asymmetric reaction fronts Ri (^) will attain nega- 
tive values in the vicinity of ^ = 0. 



2.3 Reaction fronts of type III 

It turns out that y^o(0 and Bo{^) may be nonmonotonic 
functions of ^. In this case, which we call reaction front of 
type III, Ti-iiS.) has a single maximum surrounded by two 
minimums at which TZi{^) < 0. All these properties are 
clearly seen in Fig. obtained for Dc = 10, Db = 0.1, 
bo = 0.5, and g — 2. 

This type of a reaction front can be uniquely identified 
by determining whether the maximal value of a, denoted 
Omaxj exceeds ap (or, similarly, whether 6max > 

). We 
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Fig. 6. The maximal relative increase of the concentration 
of species B, Zifemax = (&max — bo) /bo, as a function of Dc for 
g = 0.01, bo = 0.01 (solid lines), 0.001 (dashed lines), Db = 
(o), 0.5 (o), and 2 (x). Arbitrary units. 

employed this criterion in our numerical calculations. On 
extensive scanning of the 4-dimensional space of free pa- 
rameters we came to the conclusion that the necessary 
and sufficient condition for this type of the asymptotic 
reaction front reads 

Dc >max{DA,DB). (17) 

However, only for Dc ^ ma.x{DA, Db) is the effect re- 
ally significant. We also found that the maximal relative 
increase in concentrations of species B, Zi6max = (^max — 
bo) /bo, is an increasing function of Dc and a decreasing 
function of both bo and Db- In particular Z\6,„ax turns 
out to be very sensitive to changes of bo, i.e. a parame- 
ter that can be easily controlled experimentally. As for g, 
Z\6niax attains a maximal value at (7 « 1 and decreases as 
g —> or g 00. These findings are depicted in Fig. IHl 
which presents Zi6max as a function of Dc obtained for 
g = 0.01, bo = 0.01 (solid lines), 0.001 (dashed lines), and 
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Db = (circles), 0.5 (diamonds), 2.0 (crosses). As we can 
see, /i6max can attain quite high values, exceeding 60%. 

The unusual features of this asymptotic solution can 
be explained as follows. For Dc S> Da-Db molecules C 
quickly diffuse away from the reaction layer. They may 
thus form a region where the reverse reaction dominates 
the forward one, leading to 7?.i(^) < 0. The same phe- 
nomenon brings about production of additional backward 
reaction products A and B outside the main reaction zone. 
For suitably chosen system parameters this can result in a 
situation where Ao{£,) and Bo{S.) are nonmonotonic. This 
effect should become more pronounced with increased ve- 
locity of the reaction zone (i.e., when bo, Db — > 0) and 
becomes negligibly small as g (negligible backward 
reaction) or 5 — > oo (negligible forward reaction). 



3 Conclusions 

We have analyzed numerically the large-time properties 
of reaction fronts formed in reversible reaction-diffusion 
systems of type A -I- B ^ C. We found that, depending 
on the values of control parameters, reversible reaction 
fronts can be divided into three categories. In reaction 
fronts of type I the local reaction rate is always positive 
and has a well defined, single maximum. In reaction fronts 
of type II the local reaction rate has two maxima, moving 
in opposite directions, and a single minimum, which can 
attain a negative value. In reaction fronts of type III the 
local reaction rate has a single maximum surrounded by 
two minima, at which it attains negative values; moreover, 
the concentrations of species A and B are here locally 
larger than their initial values. Our numerical calculations 
indicate that the condition for this type of reaction front 
is given by a formula Dc > nia.x{DA, Db) and that the 
effect of the local increase in concentration of species A or 
B can be easily controlled experimentally through their 
initial concentrations oq or 69. 

We found that the large-time behaviour of reversible 
reaction-diffusion systems is richer than that of irreversible 
ones. Depending on the values of the control parame- 
ters one can expect qualitatively different asymptotic solu- 
tions. Although the "anomalous" effects are rather small, 
we believe that they could be observed experimentally. It 
would be particularly interesting to investigate effects of 
nonmonotonic dependence of concentrations of reactants 
A and B on the space coordinate x in systems with reac- 
tion fronts of type III. If the A -I- B ^ C reaction were a 
part of a more complex reaction scheme such that addi- 
tional reaction steps (e.g. precipitation) could occur only 
above some threshold values of the reactant concentra- 
tions (e.g. nucleation thresholds), setting oq or just be- 
low such a threshold value might lead to some interesting 
phenomena. An example of such a complex process is for- 
mation of the Liesegang patterns, which are quasiperiodic 
precipitation patterns emerging in the wake of a mobile 
chemical reaction front ^24..,25j . Our study indicates that 
it should be possible to obtain similar precipitation pat- 
terns of species B in reaction-diffusion systems with re- 



versible reaction of type A -f B ^ C, diffusion coefhcients 
Dc S> Da, Db, and initial concentrations uq ^ bo. 
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